function[y] = F_diff(p_star)
global gamma p1 rho1 c1 p2 rho2 c2 u1 u2;

if(p_star >= p1)
    f1_diff = (( (gamma+1)/2 *(p_star/p1) + (gamma-1)/2/gamma )^0.5 - ...
        (p_star-p1)/2*( (gamma+1)/2/gamma *(p_star/p1) + (gamma-1)/2/gamma)^(-0.5) * (gamma+1)/2/gamma/p1)...
        /(rho1*c1*((gamma+1)/2/gamma *p_star/p1  + (gamma-1)/2/gamma));
end
if(p_star <p1)
    f1_diff = c1/gamma*(p_star/p1)^((gamma-1)/2/gamma)*(p_star)^(-1);
end
if(p_star>=p2)
    f2_diff = (( (gamma+1)/2 *(p_star/p2) + (gamma-1)/2/gamma )^0.5 - ...
        (p_star-p2)/2*( (gamma+1)/2/gamma *(p_star/p2) + (gamma-1)/2/gamma)^(-0.5) * (gamma+1)/2/gamma/p2)...
        /(rho2*c2*((gamma+1)/2/gamma *p_star/p2  + (gamma-1)/2/gamma));
end
if(p_star <p2)
    f2_diff = c2/gamma*(p_star/p2)^((gamma-1)/2/gamma)*(p_star)^(-1);
end
y = f1_diff+f2_diff;
end